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ABSTRACT 

A recent analysis of data collected by the Planck satellite detected a net dipole at the location of X-ray selected galaxy 
clusters, corresponding to a large-scale bulk flow extending at least to z ~ 0.18, the median redshift of the cluster 
sample. The amplitude of this flow, as measured with Planck, is consistent with earlier findings based on data from 
the Wilkinson Microwave Anisotropy Probe (WMAP). However, the uncertainty assigned to the dipole by the Planck 
team is much larger than that found in the WMAP studies, leading the authors of the Planck study to conclude that 
the observed bulk flow is not statistical significant. 

We here show that two of the three implementations of random sampling used in the error analysis of the Planck study 
lead to systematic overestimates in the uncertainty of the measured dipole. The first method, rotation around the 
Galactic pole (the Z axis), increases the uncertainty of the X and Y components of the dipole and artificially reduces 
the significance of the dipole detection from 98-99% to less than 90% confidence. The second method, utilizing random 
simulations of the sky, does not take into account that the actual realization of the sky leads to filtered data that have 
a 12% lower root-mean-square dispersion than the average simulation. When either effect is taken into account, the 
corrected errors agree with those obtained with the third estimation method, and the resulting statistical significance 
of the dipole measured by Planck is consistent with that of the WMAP results. 
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The intrinsic difficulty of determining peculiar velocities 
from galaxy redshifts and distance indicators led Kash- 
linsky & Atrio-Barandela (2000) to propose an alternative 
method of probing the velocity field on large scales. Galaxy 
clusters leave an imprint on the Cosmic Microwave Back- 
ground (CMB) in the form of distortions in the CMB black- 
body spectrum caused by the Sunyaev-Zel'dovich (SZ) ef- 
fect (Sunyaev & Zel'dovich 1970, 1972). Two different 
mechanisms contribute to the SZ effect: the thermal com- 
ponent (TSZ) is caused by the thermal motions of electrons 
in the potential wells of clusters, whereas the kinematic 
component (KSZ) is due to the motion of the cluster as a 
whole. We noted that any bulk flow of clusters would pro- 
duce a dipole in the anisotropy temperature in the direc- 
tion of clusters. Since this signal is small compared to the 
sampling variance of the intrinsic CMB dipole at the same 
positions, we proposed to use the statistical properties of 
the CMB data to filter out the dominant cosmological com- 
ponent, thereby increasing the signal-to-noise ratio of any 
contribution from a bulk-flow dipole. 

In Atrio-Barandela et al. (2008) and Kashlinsky et al. 
(2008, 2009) we presented results of our application of 
this method to Wilkinson Microwave Anisotropy Probe 
(WMAP) 3-year data, later extended to WMAP 5-year 
and 7-year data (Kashlinsky et al. 2010; Kashlinsky, Atrio- 
Barandela & Ebeling 2011). For a sample of ^700 X-ray- 
selected clusters, we detected a persistent dipole, measured 
at the cluster positions within apertures (of 25' radius) that 



contain zero TSZ monopole. The dipole is roughly aligned 
with the CMB dipole and can be traced to cluster redshifts 
exceeding z ~ 0.2; its amplitude correlates with that of the 
monopole within apertures of 10' radius. We interpreted 
this signal, associated exclusively with clusters, as evidence 
of a large-scale flow of amplitude ~ 600-1000 km s _1 that 
might encompass the whole observable horizon. To illus- 
trate the significance of this result we note that a "Dark 
Flow" of this amplitude, if real, would be equivalent to the 
all-sky CMB dipole being primarily of primordial origin, 
intrinsic to last scattering surface. 

Our theoretical and numerical estimates indicated that 
our dipole detection is significant at the 99.4% confidence 
level. However, independent confirmation of this result is 
still lacking, and several studies have challenged our re- 
sults. Keisler (2009) confirmed the existence of the dipole 
detected by us, but claimed that it is not statistically signif- 
icant. It was shown though (Atrio-Barandela et al. 2010) 
that Keisler neglected to subtract the dipole outside the 
Galactic mask, causing the error bars of his measurement 
to be overestimated. More recently, Osborne et al. (2011) 
and Mody & Hajian (2012) found no bulk flows in WMAP 
data in analyses based on filtering schemes that differ from 
ours. However, both teams of authors implicitly assumed 
that clusters have the same angular extent in the origi- 
nal as in the filtered data, and that the electron pressure 
and the electron density in clusters follow the same radial 
profile. Both assumptions are incorrect and render these 
filters insensitive to bulk flows. In Atrio-Barandela et al. 
(2013) we demonstrated that correct implementation of ei- 
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ther of these alternative filtering schemes leads to results 
that are consistent with ours. Two recent studies using 
galaxies, rather than clusters, to search for the KSZ sig- 
nature of a large-scale bulk flow in CMB data have also 
challenged our findings. The CMB-galaxy cross-correlation 
study by Li, Zhang & Chen (2012) found a statistically 
anisotropic component in the CMB on scales out to z ~ 0.2, 
but ruled out the Dark Flow as its cause, due to the very 
high galaxy velocities implied by this interpretation. We 
note though that peculiar velocities in perfect agreement 
with our findings would result if their cross-correlation sig- 
nal were dominated by cluster galaxies rather than field 
galaxies. Lavaux, Afshordi & Hudson (2013) aimed to use 
the KSZ effect from the plasma halo of galaxies to probe 
the peculiar velocity field and found results on small scales 
(<50/i _1 Mpc) that are consistent with those obtained us- 
ing galaxy distance indicators. On scales of ^500/i _1 Mpc, 
however, their analysis leads to an upper limit of 470 km 
s _1 (95% significance) for the bulk-flow velocity, in conflict 
with the Dark Flow amplitude measured by us. It should 
be noted though that these conclusions depend sensitively 
on the spatial distribution of ionized gas in and around 
galaxies, which is presently ill-constrained by observations. 

A conclusive determination of the Dark Flow is expected 
to emerge from an authoritative study based on CMB data 
collected by the Planck satellite. With its higher resolution, 
better frequency coverage, and lower noise levels Planck 
is much better suited for studies of the SZ effect than 
WMAP. A comparison between the dipole measurements 
from the two missions will also allow a much-needed con- 
sistency check, since differences between the instruments 
and scanning strategies result in data with different sys- 
tematics. In this context, the Planck Collaboration has 
recently reported (Planck Collaboration. Planck 2013 re- 
sults. XXVII. 2013) the aberration and modulation effects 
due to the motion of the Sun with respect to the CMB 
frame and found a significant agreement with the Doppler 
boosting of the CMB monopole. However, this measure- 
ment does not constrain the motion of the Local Group 
but the amplitude of its component in the direction of the 
solar motion. They have also reported the first Planck- 
based KSZ measurement of the peculiar motions of clusters 
(Planck Collaboration 2013, Planck Intermediate Results 
XIII, hereafter PIR-13) using their own catalog of X-ray 
selected clusters (Piffaretti et al. 2011) and concluded that 
the Planck constraints on peculiar velocities were consistent 
with the concordance ACDM model. 

In this paper we analyze and evaluate the error estima- 
tors employed by PIR-13. We begin by summarizing the 
PIR-13 results and methods in Sec. 2. Sec. 3 briefly re- 
views the formalism regarding error estimation developed 
by Atrio-Barandela et al. (2010), before we discuss in detail 
two problematic implementations of error estimation used 
in PIR-13, rotations and numerical simulations, in Sec. 4 
and 5, respectively. We present our conclusions in Sec. 6. 

2. Results and methods of the first Planck 
bulk-flow study 

The PIR-13 results characterize the average peculiar veloc- 
ities of clusters at different depths using three complemen- 
tary diagnostics: the radial peculiar-velocity average, the 
root-mean-square (rms) velocity, and the bulk-flow velocity. 
The underlying analysis employs three different filtering 



schemes, including the one introduced by us (Kashlinsky 
et al. 2009). The study's findings for all three diagnostics 
mentioned above are consistent with the very low ACDM 
predictions for the average peculiar velocities of galaxy clus- 
ters. Specifically, the PIR-13 abstract asserts that "there is 
no detection of [a] bulk flow as measured in any comoving 
sphere extending to the maximum redshift covered by the 
cluster sample". To fully appreciate the nature and signif- 
icance of this assertion, it is critical to understand that, 
when using our filtering scheme, PIR-13 detects the same 
dipole signal at cluster positions as previously found by 
us in WMAP data (PIR-13, Fig 10). What renders their 
bulk-flow detection insignificant is thus not the absence of 
a dipole signal, but the uncertainty assigned to it. Before 
analyzing PIR-13's error estimates in more detail in the fol- 
lowing section, we here briefly discuss several key aspects 
of their work. 

The PIR-13 study is based on the Meta Catalogue of X- 
ray detected Clusters of galaxies (MCXC; Piffaretti et al. 
2011), a compilation of publicly available X-ray cluster sam- 
ples drawn primarily from ROSAT All-Sky Survey projects, 
but also including mostly low-mass systems from serendip- 
itous X-ray cluster surveys. Comprising 1750 clusters and 
rich galaxy groups, the MCXC is comparable in size to our 
proprietary catalog of 1558 X-ray selected clusters, but con- 
tains more low-luminosity systems and fewer very X-ray lu- 
minous clusters. By excluding MCXC entries close to bright 
point sources in any single-frequency Planck map, clusters 
lying in regions of high Galactic emission, and systems with 
estimated masses below 1O 13 M , the sample used by PIR- 
13 is reduced to 1405 clusters. PIR-13 conducted their anal- 
ysis using a CMB map that is the weighted average of the 
single-frequency Planck maps. This so-called 2D Internal 
Linear Combination (2D-ILC) map introduces additional 
masking, resulting in a final sample of 1321 MCXC clusters. 
For comparison, our sample was limited to the 1205 clusters 
with X-ray luminosities L x [0.l - 2.4 keV] > 0.5 x 10 44 erg 
s _1 , and our results were obtained from 796 clusters that 
fall outside the extended KQ75 mask. 

Comparing the WMAP and Planck cluster dipoles al- 
lows us to identify and account for previously undiagnosed 
systematic effects. Because of such systematic differences it 
is possible that, even if the measured values are the same, 
the error bars derived from WMAP and Planck data differ. 
The main reasons for such differences are twofold. On the 
one hand, the measurement uncertainties are dominated by 
the residual cosmological anisotropy left by the filter, com- 
mon to all frequencies. Due to its higher angular resolution, 
Planck probes the CMB spectrum up to i ~ 2000, com- 
pared to I 1000 for our WMAP filtered maps. Planck's 
filtered maps will thus contain a larger CMB residual than 
those based on WMAP data, which could lead to larger er- 
rors and lower significance of the KSZ dipole detection. On 
the other hand, a larger cluster catalog must improve the 
significance of the result (assuming that the cluster mass 
range is comparable). Further systematic differences be- 
tween Planck- and WMAP-based analyses of this kind in- 
clude that the TSZ contribution can be determined and 
monitored much better with Planck data by comparing 
measurements at frequencies above and below the TSZ null 
at 217 GHz. It follows that any correlation between the 
TSZ monopole within 10' (radius) apertures and the KSZ 
dipole within 25' apertures (Kashlinsky et al. 2010, Kashlin- 
sky et al. 2011) can also be established better with Planck 
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than with WMAP. These obvious tests present some dif- 
ficulties though: since, in our WMAP-based analysis, the 
dipole signal was measured within a radius of 25' around 
each cluster, the high-frequency component of the CMB 
residual probed only by Planck's higher angular resolution 
would be erased. In addition, the correlation between the 
SZ monopole and dipole observed by us indicates that low- 
mass systems contribute little to the signal. A larger frac- 
tion of low-mass clusters in the MCXC could thus dilute 
the dipole, although we believe those are not the reasons 
for which PIR-13 failed to detect a significant cluster dipole, 
as we shall show below. 

PIR-13 applied several filtering schemes in Legendre 
space. While the actual filters differ from those previously 
used and described in the literature, the PIR-13 filter im- 
plementations share the drawbacks of the filters used by 
Osborne et al. (2011) and Mody & Hajian (2012). Specifi- 
cally, the PIR-13 study assumes that clusters have the same 
size in the real as in the filtered data, and adopts isothermal 
cluster models to calibrate these filters. Both of these as- 
sumptions, however, render filters insensitive to bulk flows, 
as shown by Atrio-Barandela et al. (2013). In addition, 
PIR-13 uses the Aperture Photometry filter, first intro- 
duced in the context of bulk-flow measurements by Kash- 
linsky & Atrio-Barandela (2000). Contrary to other filters, 
the Aperture Photometry filter operates in real space: the 
intrinsic CMB signal is removed by subtracting the temper- 
ature averaged within an annulus (defined by radii 9i n , O ut) 
from the temperature averaged within the inner disc of ra- 
dius Bi n . While this process indeed removes the CMB, it 
also removes part of the SZ anisotropy of the cluster. In 
order to estimate the missing fraction, the profiles of the 
electron pressure and electron density need to be known 
for each cluster. Assuming isothermality simplifies the task, 
but biases the result. If the electron density falls less rapidly 
than the electron pressure, then the KSZ component is re- 
moved more efficiently than the TSZ component, reducing 
the cluster dipole. In addition, all cluster annuli must have 
the same radii if the CMB is to be removed uniformly across 
the sky, but PIR-13 defines different radii for each cluster, 
causing the residual CMB to vary from cluster to cluster. 
A detailed discussion of the resulting complications will be 
given in a forthcoming paper. 

Finally, and most importantly, the PIR-13 analysis uses 
three methods to estimate the uncertainties of their mea- 
surements: (1) dipole measurements performed on the fil- 
tered maps at randomized cluster positions, (2) dipole mea- 
surements performed on random realizations of the filtered 
data, but at the real cluster positions, and (3) dipole mea- 
surements performed by rotating the template of cluster 
positions on the real data. PIR-13 offers no discussion of 
these estimators' relative efficiency or biases, which is sur- 
prising in view of the fact that the resulting error estimates 
differ greatly. When using (1), the significance of the KSZ 
dipole found by PIR-13 using Planck data is higher than 
the one found by us from WMAP data, but it is lower when 
using methods (2) and (3). A bias inherent in method (2) 
is well established and due to the power spectrum of the 
residual CMB component in the actual realization of the 
sky being about 10% smaller than that of the average ran- 
dom realization of the sky (Atrio-Barandela et al. 2010). In 
the following sections, we discuss key aspects of the error 
computation as well as the shortcomings of PIR-13's error 
estimation methods (2) and (3) in more detail. 



3. On the Error Bars of the Bulk Flow 
Measurement 

Since clusters subtend a small solid angle on the sky, sam- 
pling variance is the largest source of uncertainty in the 
determination of bulk flows (Kashlinksy & Atrio-Barandela 
2000). Random dipoles from the cosmological CMB com- 
ponent have much larger amplitude than any KSZ dipole 
caused by the bulk motion of a cluster sample. To sup- 
press this noise term, we suggested to filter out the CMB 
signal using the statistical properties of the temperature 
anisotropy field. In addition, contributions from the Galaxy 
and point sources are removed by application of suitably de- 
fined masks. We assume use of the extended WMAP 7yr 
mask (the KQ75 mask) in the following description of the 
convolution process that creates the filtered maps: 

1. Each CMB map is multiplied by the KQ75 mask. 
Monopole, dipole, and quadrupole are removed from the 
remaining pixels. 

2. The resulting CMB map is expanded into spherical har- 
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). The radiation power 
//sky, where the division 



by /sky, the fraction of the sky outside the KQ75 mask, 
accounts for the power removed by the masking process. 

3. The power spectrum is filtered, and a Legendre trans- 
formation back to real space is applied: AT m = 
J2 a em FeYim- We used the following Wiener- type fil- 
ter: F t = [Cf y - Cf CMB B 2 £ ]/Cf y , where B t is the 
antenna beam, and C^ CMB is the theoretical CMB ra- 
diation power spectrum of the concordance model that 
best fits the data. 

4. Finally, the filtered map AT m is multiplied by the mask 
and any monopole and dipole introduced by the filtering 
process are removed from the remaining pixels. 

In Atrio-Barandela et al (2010) we presented a compre- 
hensive discussion of the errors associated with our method 
and showed that the variance of the monopole and dipole 
components is 



Var(a ) 



°CMB,fil , a noise,fil 
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with (m) = (X,Y,Z) = (cos I cos 6, sin I cos b, sin b) being 
the direction cosines of the cluster positions on the sky. In 
this expression, cr n oise,fii and <7cMB,fii are the residual CMB 
and noise in the filtered maps. The noise is uncorrelated 
from pixel to pixel, whereas the correlation function of the 
residual CMB noise crosses zero at 6 ~ 1°. Therefore, the 
sample variance of the contribution from the residual CMB 
scales approximately as the number of clusters N c \, while 
the noise contribution scales with the number of pixels 7V p i X . 

When estimating measurement uncertainties, we have 
to take into account (a) that clusters are not randomly dis- 
tributed in the sky and that any anisotropy in their dis- 
tribution could increase the errors, and (b) that, in ad- 
dition to instrumental noise, filtered maps also contain un- 
known residuals of foreground emission and other systemat- 
ics that cannot be precisely modeled. To take into account 
all sources of statistical uncertainty, error bars would have 
to estimated using both the filtered data and the actual 
template of clusters in the sky. This is not possible (see 
also Sec. 2]). However, errors computed using random tem- 
plates on the filtered data showed a behaviour very close to 
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the prediction of eq. ([T]), indicating that foreground resid- 
uals and other possible systematics can be safely neglected 
(Atrio-Barandela et al. 2010; Kashlinsky, Atrio-Barandela 
& Ebeling 2012). To further explore item (a), the impact 
of anisotropics in the cluster distribution, we will here con- 
sider two templates: Tl, defined as the 796 clusters in our 
catalog with L x > 0.5 x 10 44 erg s" 1 in the ROSAT broad 
band (0.1-2.4 keV), and T2, a subset of Tl, comprising only 
the 327 most X-ray luminous clusters (Lx > 2.0 x 10 44 erg 
s" 1 ) at z < 0.3. 

Convolving the cluster profile with the filter increases 
the size of clusters in the filtered map (Atrio-Barandela et 
al. 2013). We assign all clusters the same angular radius 
of 25', which roughly corresponds to the aperture within 
which the monopole vanishes for WMAP 7-year data. At 
Healpix resolution A sidc = 512 (Gorski et al. 2005), sam- 
ples Tl and T2 occupy 32,700 and 13,500 pixels, respec- 
tively (about 1% and 0.4% of the sky). Since gravitational 
instability drives low-mass systems towards the more mas- 
sive clusters, we expect the distribution of the Tl sample 
to be less isotropic than that of the T2 template, offering 
the possibility to quantify the impact of anisotropics on the 
measurement errors. Using those templates and WMAP 7- 
year data from the Wl Differencing Assembly (DA) , we find 
fCMB.fii ~ 30/iK and cr no isc.fii ~ 75/iK. Given the relatively 
modest number of clusters in our samples Tl and T2, eq. ([1} 
then indicates that the residual CMB in the filtered map 
dominates the errors, as suggested independently by Keisler 
et al. (2009) and Kashlinsky et al. (2010). Simulations of 
filtered maps thus need to contain only the CMB, and can 
neglect noise, foreground residuals, and other systematics. 

4. Error estimates from rotated cluster templates 

As mentioned before, a rigorous computation of errors 
would require measuring dipoles at random locations in the 
filtered data, outside the cluster apertures, but using the 
same template that describes the real cluster positions used 
in the analysis. Because of the complex geometry of the 
mask, this is not possible. In Atrio-Barandela et al. (2010) 
we discussed the methods labelled (1) and (2) in Sec. [2j 
option (1) entails measuring dipoles in the filtered data, 
using random cluster templates that have no overlap with 
the real one, whereas option (2) measures dipoles in simu- 
lated data, but at the position of the real clusters. When 
applying method (1), we account for the contributions to 
the measurement error from foreground residuals and other 
systematics as far as they are known; with method (2), we 
account for the effect of anisotropic distribution of our clus- 
ter template on the sky as detailed later in this section. In 
our simulations we do not find any systematic differences 
between results obtained for real and simulated cluster tem- 
plates beyond cosmic variance. This is understandable since 
clusters, at random or real positions, are, by design, selected 
outside the Galactic plane, and hence all templates share 
the large-scale inhomogeneity of the mask. 

In PIR-13 a third method was introduced: (3) dipoles 
were measured in the filtered data at locations obtained 
by rotating the cluster template around the Z axis in steps 
of one degree. The Planck Collaboration expected this 
method to be optimal, since it includes all possible sys- 
tematics of both the filtered data and the cluster template, 
including anisotropics and angular correlations. Unfortu- 
nately, this method is very inefficient and yields system- 



atically larger errors than homogeneous random sampling. 
For one, the fixed step size of one degree allows only 359 dif- 
ferent measurements and thus yields very few independent 
dipole estimates. In addition, a rotation does not move the 
clusters within the template alike. For instance, under rota- 
tion about the Z axis, the Coma cluster, located at b = 88°, 
never moves by more than four degrees from its initial posi- 
tion. Finally, for small rotation angles, nearby clusters that 
are very extended on the sky will continue to contribute to 
the random dipole, correlating the measurements. As a re- 
sult, the space of all possible dipoles is poorly sampled; the 
measured random dipoles overpopulate the tail of the dis- 
tribution, leading to overestimated uncertainties that dilute 
the significance of any real result. 

In the following we discuss these systematic effects in 
more detail. 

4.1. Artificial correlations 

Fig. [U shows the monopole, dipole, and angular direction 
of the dipole measured by rotating the Tl template (796 
clusters) in the same manner as prescribed by PIR-13. We 
show results obtained with the W4 DA since it yields the 
largest dipole measured in the WMAP W-band, making 
it easier to identify any trends and biases. While, by de- 
sign, all clusters in the unrotated Tl template fall outside 
the KQ75 mask, this is no longer the case after rotation^. 
The ensuing loss of clusters due to rotation (up to 10% of 
all pixels fall outside the mask for a given rotation angle) 
introduces additional variance in the measured dipoles. 

Fig. Q] illustrates several of the artifacts introduced by 
PIR-13's rotation method. For instance, the measured 
monopole appears to repeat with a ~ 120° period, and 
the distribution of dipoles is similarly inhomogeneous. If a 
rotated template indeed sampled random dipoles, the angu- 
lar directions of these dipoles would have to be distributed 
randomly. As shown by Fig. [T] this is clearly not the case; 
strong correlations are visible, especially for the Galactic 
longitude I of the dipole direction. Moreover, the ampli- 
tudes of the monopole and dipole are also correlated. For 
the Tl template, the correlation matrix of the dipole com- 
ponents is: 

, x / 1 0.05 -0.13 \ 

C(au, aij )= 2 %ffi W2 = °-° 5 1 ~°- 26 

(a H ) 1 \a Xj ) i y _ 13 _ 26 i J 

(2) 

In this expression, = (X,Y,Z). Remarkably, the 

dipole component along the rotation axis, Z, correlates 
more strongly with the other two components than X and 
Y. These correlations are larger than the error bars esti- 
mated from method (1), which uses the real data and ran- 
dom cluster templates, and where for 400 dipoles the off- 
diagonal terms never exceed 7% of the diagonal terms (see 
Kashlinsky, Atrio-Barandela & Ebeling, 2012). 

From the results shown in the top right panel of Fig. [T] 
we can immediately compute the significance of the dipole 



In principle the cluster template could of course be rotated 
about any axis but, due to the symmetry of the KQ75 mask 
with respect to the plane of the Galaxy, the choice of an axis 
other than the one through the Galactic poles would excise a 
much larger fraction of pixels for a rotated cluster template. 
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measurement. We find the dipole measured for the unro- 
tated Tl template to be significant at the 92% confidence 
level; for the T2 subsample, the significance increases to 
96%. Large variations in the significances thus derived are 
expected, owing to the small number of dipole measure- 
ments allowed by this method. Indeed, PIR-13 report yet 
another value (89%) when using our filter and their MCXC 
cluster template. 

4.2. Inflation of X and Y uncertainties 

When the template is rotated, clusters do not move homo- 
geneously on the celestial sphere. For rotations about the 
Z axis, the number of clusters at constant Galactic lati- 
tude, b, is fixed. Since clusters move more slowly close to 
the Galactic pole than they do near the Galactic plane, the 
Z component of the dipole is more heavily sampled than 
the X and Y components, resulting in smaller uncertainties 
for the former than for the latter. We can quantify this 
over- and undersampling by using eq. ([T]), which gives the 
error of each dipole component in units of the error of the 
monopole. 

For clusters that are isotropically distributed on the sky 
the term (nf) in eq. (TjTJ equals 1/3. This is the minimum 
variance estimator since the three quantities an are derived 
from the same data set as hq, and the error of each dipole 
component is given by a(au) = \/^a{ao). In Table 1 we list, 
for all three error-estimation methods outlined before, the 
fractional increase of the errors caused by the anisotropy of 
both the cluster distribution and the sampling of the sky, 
relative to uniform sampling. Note that these error esti- 
mates only include the effect of the geometry of the cluster 
template, but not the additional variances introduced by 
the limited number of dipole estimates and from the sys- 
tematics intrinsic to the data. Hence, the listed Figs, are 
the minimum error that can be achieved for a given cluster 
configuration. 

Note that, even for perfectly random sampling, i.e., 
method (1), the presence of the mask causes the error in 
the X and Y components to increase with respect to homo- 
geneous sampling of the full celestial sphere, while the error 
in the Z component decreases. The same effect is evident 
for method (2) which results in slightly larger uncertainties 
due to the increased variance introduced by the anisotropy 
of the distribution of clusters in the sky. The by far largest 
deviations from uniform sampling, however, are induced by 
method (3) which employs rotations of the cluster template. 
Note (rightmost two columns in Table 1) that method (3) 
performs worst in the Y direction where it leads to a 15% 
increase in the uncertainty compared to method (2). It is 
the Y component of the dipole that dominates the Dark 
Flow signal (Kashlinsky et al. 2010). 

Fig. [2] shows a stack of the Tl template and its 359 
one-degree rotations, illustrating the non-uniform sampling 
resulting from PIR-13's rotation method. The observed 
banding, predominantly at high Galactic latitude, is a con- 
sequence of the differences in both the number of clusters 
for different values of b and in the angular displacement of 
clusters in Galactic longitude as a result of the rotation. 
The celestial sphere is sampled inhomogeneously, increas- 
ing the errors and resulting in dipoles that are strongly 
correlated (see also eq. ©, and Fig. |T|) . 



The uncertainties resulting from the three different er- 
ror estimators (Table 1) can be translated into significances 
of detection for a dipole signal of given amplitude. Fig. [3] 
shows these significances, based on 359 realizations for each 
of the three methods. A dipole of amplitude 3.8/xK would 
be significant at the 97% confidence level if the error es- 
timate were obtained from 359 measurements of random 
dipoles that sample the celestial sphere perfectly uniformly. 
For methods (1), (2), and (3), applied outside the KQ75 
mask, this number is 96%, 94% and less than 90%. In re- 
ality, the loss of significance caused by method (3) is even 
higher, due to the correlation between the dipole compo- 
nents, not accounted for in these theoretical estimates. 

In summary, method (3) (sampling of the real CMB sky 
by rotating the cluster template about the Z axis) clearly 
and systematically degrades the significance of a dipole 
measurement compared to methods (1) and (2). Its sam- 
pling of the celestial sphere is the most inefficient and most 
anisotropic of the three estimators, and the method is es- 
pecially insensitive to the Y-component of the dipole. 

5. Error estimates from CMB simulations 

Although method (3), rotation, is the most inefficient and 
biased error estimation technique discussed here, care also 
needs to be taken when interpreting the errors obtained 
from numerical simulations of the CMB sky, i.e., method 
(2). In general, simulated maps feature larger variances 
than the real sky and consequently lead to larger error es- 
timates which, in turn, cause the significance of any real 
signal to be underestimated. We have identified three ef- 
fects that contribute to this bias: (a) in the filtered data the 
power leaks from regions outside the mask and reduces the 
variance on pixels outside the mask compared with a simu- 
lated map (b) the power spectrum of the actual realization 
of the sky C| ky F^ 2 is smaller than the power spectrum of the 
theoretical ACDM concordance model Cf CT>M BfFf used 
to generate the simulated sky, and finally, (c) accidental 
sampling of a region with non-zero dipole, if the monopole 
and dipole were not removed from the region outside the 
KQ75 mask, prior to computing the dipole for a random 
sky. 

To quantify the impact of power leakage to the masked- 
out part of the sky we generated 1,000 realizations of the 
CMB sky corresponding to the specifics of WMAP 7- year 
data Wl DA. Simulated maps were computed using Healpix 
N s ide = 512 resolution, and multipoles were sampled up to 
(max — 1024. The simulated (a^ m ) coefficients were multi- 
plied by the Wl beam before transformation to real space. 
Next, each map was multiplied by the mask, converted back 
into spherical harmonics, filtered, and transformed once 
more into real space as described in Sec. [31 For each real- 
ization, the Wiener-type filter used has the same functional 
form as specified in Sec. [3J but C| ky is now replaced by the 
power spectrum of the simulated "data", cf m ~ sky + Cf olse . 

Here Cf m ' sy is the power spectrum of each simulated sky, 
corrected, if necessary, for the fraction of power removed 
by the mask, and Cf olsc is the power spectrum of the noise. 
For each simulation we computed the rms of the filtered 
map inside and outside the mask. If in the initial data we 
do not multiply by the mask, the filtered map has an aver- 
age power of a — 30±3/iK. If in the initial data we multiply 
by the KQ75 mask, then the power on the fraction of the 
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Table 1. Relative uncertainty in the measurement of the three spatial components of the dipole, obtained with different error- 
estimation methods. Values of unity are obtained for a uniform sampling of the sky. The error-estimation methods are labeled as 
before: (1) real data, random template; (2) simulated data, cluster template; (3) real data, rotated cluster template. All methods 
consider only positions outside the mask. 





Method (1) 


Method (2) 
Tl T2 


Method (3) 
Rl R2 


Ratio 
Rl/Tl R2/T2 


[{JP)/3]-V* 
[{Y 2 )/3}-^ 2 
[<^)/3]-Va 


1.124 
1.060 
0.871 


1.207 1.192 
1.083 1.093 
0.883 0.828 


1.279 1.302 
1.235 1.259 
0.754 0.749 


1.060 1.092 
1.140 1.152 
0.912 0.906 



mask outside the mask was a out = 28 ± 2/iK, while in the 
region excluded by the mask it was a m — 15 ± 3/j.K. This 
effect of power leakage is illustrated in Fig. Q] for one of our 
simulations. We found the average fraction of power leaked 
to the mask to be 0.061 ± 0.016. Therefore, use of filtered 
maps generated using realizations of the filtered power spec- 
tra not taking onto account the effect of the mask will, on 
average, increase the errors by ~ 6%. 

The second effect contributing to the artificially high 
variance of the simulated CMB sky is illustrated in Fig. [5] 
where we plot the mean and rms dispersion of the filtered 
power spectrum of the simulations. We define the cumula- 
tive variance out to a given £ as 

1 1 

<K0 2 = ^£(2i+i)c< ( 3 ) 

where Cj is the filtered power spectrum. In Fig. [5] the solid 
black line represents the filtered power spectrum from the 
Wl DA data that contains both residual CMB and noise 
(solid black line), the thin red line represents (^ ACDM ' fil = 
Cf c ™B}{Ff 1 ) 2 , where C^ ACDM is the theoretical power 
spectrum of the ACDM cosmology and F^ 1 the filter con- 
structed with the Wl data. Also shown is the filtered 
power spectrum, averaged over 1,000 simulations, (dot- 
dashed blue line) and the rms dispersion around this mean 
(dashed blue lines). Let us remark that a(£ ma x) is slightly 
larger than crcMB,fii of eq. JT]). We assign all clusters a ra- 
dial extent of 25', spatial scale that corresponds to £ ~ 400. 
Then, CMB residuals beyond £ ~ 400 do not contribute to 
the variance of the monopole and dipole measured within 
cluster apertures. As demonstrated in Atrio-Barandela et 
al. (2010), (j cMB,fii is approximately 15— 17/LtK, the ampli- 
tude of a(£) for 't = 300 - 400. Note that while in Fig. 
£,ACDM,fii j g yer y s j m ji ar to the average power spectrum 
of one thousand simulated maps (filtered using the same 
pipeline as applied to the data), the power of the actual 
realization of the sky Cg ky (F^ 1 ) 2 is about 10% smaller at 
£ = 300 (at larger £'s, it is larger due to the noise con- 
tribution). We thus reach the same conclusion as Atrio- 
Barandela et al. (2010); adding in quadrature the power 
leak due to the mask, simulations of filtered maps over- 
estimate the errors on the monopole and the three dipole 
components by 12% with respect to the actual realization 
of the sky. 

The third effect to consider when using simulated maps 
for the estimation of errors in the dipole measurement has 
already been mentioned in Sec. [31 For each of our thousand 
filtered maps, the dipole outside the KQ75 mask must be re- 
moved before computing the dipole at the cluster locations, 
just like for the real data. This step is essential to ensure 
that the random dipoles are computed for a zero-dipole sky. 



Keisler (2009) overlooked this crucial step, thereby artifi- 
cially increasing the error bars. In Atrio-Barandela et al. 
(2010) we quantified that this omission leads to an increase 
in the error of the monopole close to 30%, of 15 — 20% for 
the X and Y components of the dipole, and of 2% for the 
Z component of the dipole. 

To quantify the impact of the aforementioned biases 
on measurements of the statistical significance of the Dark 
Flow, we applied all three error estimation methods to Wl 
filtered data. To test the relevance of anisotropics in the 
cluster distribution, we conducted the analysis for the pre- 
viously introduced Tl and T2 templates, but, for compar- 
ison, also for two randomly generated cluster templates SI 
and S2 which comprise the same number of clusters as Tl 
and T2. We generated 10 simulated maps of the sky and, 
for each of them, computed the dipole for 100 random po- 
sitions of each template, obtaining a total of 1,000 dipoles. 
Before computing the dipole at the cluster locations, we re- 
moved the dipole outside the KQ75 mask. For simplicity, 
we did not include any noise or foreground residuals as they 
contribute negligibly to the final error budget. The initial 
power spectrum was C ACDM ' fil = C^ GDM B 2 (F^ 1 ) 2 where 
(jACBM j g ACDM radiation power spectrum that best 
fits WMAP 7-year data, and F^ 1 is the Kashlinsky et al. 
(2008) filter for the Wl band. To test method (3) we also 
computed 360 dipoles using one single map, by rotating the 
Tl and T2 templates in one-degree steps about the Z axis. 
We kept the KQ75 mask fixed during the rotation, so that 
pixels of the cluster templates were excised by the mask in 
the same way as the data. 

The results of this exercise are shown in Fig. [5] for the 
Tl and SI cluster templates. Each method assigns different 
statistical significance to the dipole measured by Kashlin- 
sky et al. (2010), namely 98.6%, 96.9% and 85.8% confi- 
dence level for method (1), using the random template SI, 
method (2), using the real cluster template Tl, and method 
(3), using rotations of template Tl. Once we correct our 
simulations to remove the excess of power in the simulated 
sky compared with the actual realization of the sky, as 
shown in Fig. [SJ the significance increases to 99.6%, 98.4% 
and 90%, respectively, showing that the rotation method 
dilutes the significance. 



6. Conclusions 

Using Planck data and the filtering scheme of Kashlinsky 
et al. (2010), the PIR-13 team has detected the same dipole 
signal reported by Kashlinsky et al. (2010) and Kashlinsky 
et al. (2011) using WMAP data, but reaches a very different 
conclusion. Based on error estimates derived from simula- 
tions and rotations of the cluster template, the Planck Col- 
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laboration assigns the Dark Flow measurement a statistical 
significance below 90%. 

In this paper we have analysed the error-estimation 
methods used in PIR-13 and found biases that result in 
systematic overestimates of the uncertainty of the dipole 
measurement, and hence in a systematic underestimate of 
its significance. In particular the method of rotating the 
cluster template about the Z axis (devised by PIR-13 to 
account for angular correlations between clusters on the 
sky and all systematics present in the data) results in arti- 
ficially inflated errors, due to several systematic flaws of the 
method: (1) since only 359 dipole measurements are possi- 
ble, the resulting error estimate is not statistically robust; 
(2) the dipole components are strongly correlated with each 
other (see Fig [T] and eq. @; and (3) the sampling of the 
sky is highly anisotropic (see Fig[l|. Effects (2) and (3) 
result in errors for the X and Y components that are 5-15% 
larger than those obtained by the error-estimation methods 
discussed in Atrio-Barandela et al. (2010) and Kashlinksy 
et al. (2012) while decreasing the uncertainty of the Z com- 
ponent of the dipole. In addition, we have shown that nu- 
merical simulations, as implemented by PIR-13, also over- 
estimate the uncertainty of the dipole measurement. Due 
to cosmic variance, the actual realization of the sky has 
less power at t = 10 — 300 than the average of an ensem- 
ble of simulated skies (see Fig. [5]) which leads to errors 
derived from simulated CMB maps that are too high by 
about 12%. When all biases are corrected for, the signifi- 
cance at the 90% confidence level reported by PIR-13 for 
the Dark Flow increases to ~ 99% in agreement with our 
previous estimates (Atrio-Barandela et al. 2010). 

The strongest evidence that the measured dipole is gen- 
uinely associated with clusters was that it correlates with 
the TSZ signal measured at the cluster center. It is puz- 
zling that this correlation was not observed by the Planck 
Collaboration. If the correlation is not there implies that 
either WMAP and Planck data are systematically different 
at the position of massive clusters or that the flow con- 
verges and does not reach as deep as the scale probed by 
the most massive clusters. Either possibility would be a 
very important result. 
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methods. For a given dipole amplitude, method (1) yields the 
smallest errors and thus the highest significance (red dot-dashed 
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tations of the actual cluster template, creates the by far largest 
uncertainties and thus assigns the lowest significance. Both 
method (1) and (2) were applied to 359 realizations of a ran- 
dom cluster template or a simulated CMB sky, respectively, in 
order to match the fixed statistics of 359 rotations provided by 
method (3). 
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Fig. 1. Monopole and dipole components evaluated at locations obtained by rotating the Tl cluster template about the Z 
axis. The dot-dashed (blue) line shown in the upper right panel marks the amplitude measured at zero rotation. The bottom two 
panels show the direction of the dipole measured for a given rotation angle. Note the pronounced inhomogeneities in the resulting 
sampling of all possible random dipoles, as well as the 120° periodicity in the amplitude of the monopole. 



¥1 DA; POWER LEAKED TO THE MASKED PIXELS 




Fig. 4. Residual power in a simulated map. We use the 
WMAP 7-year KQ75 mask. For clarity, the power outside 
the mask is set to zero and the plot is limited to the range 
[—20, 20]/iK. Points above/below the specified range were given 
the value of the upper / lower bound. Note the power leaked to 
the masked-out region of the sky. 
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Fig. 5. Contribution to the rms dispersion of all multipoles up 
to t (see eq© for the WMAP Wl band. The thick solid (black) 
line corresponds to WMAP data, the thin solid (red) line to the 
best-fit concordance model (j^ CUM ' m t t ne dot-dashed (blue) line 
is the mean of 1,000 simulations containing only CMB, and the 
symmetric dashed (blue) lines mark the la dispersion around 
the mean. 
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Fig. 6. (a) Histograms of 1,000 dipoles measured from ten simulated filtered maps for the Wl DA. For each map, dipoles were 
computed placing the template at 100 random positions. The solid black and blue lines show the results obtained for the template 
of real clusters, Tl, and for the simulated template of randomly positioned clusters, SI, respectively. For comparison, the solid red 
line shows the distribution of 360 dipoles computed by rotating the Tl template on a simulated map. (b) As (a) but focusing on 
the 3 — 5.5/ik range and showing the Rl renormalized to a total of 1,000 to facilitate the comparison of the three error estimators. 
Panels (c) and (d) show the distribution of the dipole directions in Galactic coordinates. In all panels, the vertical dot-dashed line 
represents the measured dipole according to Kashlinsky et al. (2010). 
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